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ABSTRACT 

The dark energy dominated warm dark matter (WDM) model is a promising alterna- 
tive cosmological scenario. We explore large-scale structure formation in this paradigm. 
We do this in two different ways: with the halo model approach and with the help of an 
ensemble of high resolution ./V-body simulations. Combining these quasi-independent 
approaches, leads to a physical understanding of the important processes which shape 
the formation of structures. We take a detailed look at the halo mass function, the 
concentrations and the linear halo bias of WDM. In all cases we find interesting de- 
viations with respect to CDM. In particular, the concentration-mass relation displays 
a turnover for group scale dark matter haloes, for the case of WDM particles with 
masses of the order biwdm ~ 0.25keV. This may be interpreted as a hint for top-down 
structure formation on small scales. We implement our results into the halo model and 
find much better agreement with simulations. On small scales the WDM halo model 
now performs as well as its CDM counterpart. 
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1 INTRODUCTION 

Over the last decade the vacuum energy dominated cold 
dark matter (hereafter ACDM) scenario, has emerged as a 
standard model for cosmology. This owes largely to the com- 
bination of information from galaxy clustering surveys such 
as the 2dFGRS and SDSS with WMAP measurements of 
the temperature anisotropies in the microwave background 
( |Cole et al.|2005||Tegmark et al.|26"06[|Komatsu et al.|2011[ ). 
However, the nature of the two dark components in the 
ACDM model are still completely unknown and it is there- 
fore important to keep exploring alternative models and test 
their compatibility with observations. 

In the ACDM model the dark matter is assumed to be 
composed of heavy, cold thermal relic particles that decou- 
pled from normal matter very early in the history of the Uni- 
verse jPeeb!es|l^ 



1990 Jungman et al.||1996[ ). Whilst there is a large body of 
indirect astrophysical evidence that strongly supports CDM, 
there are some hints that it has shortcomings. Firstly, CDM 
galaxy haloes contain a huge number of subhaloes ( |Moore| 
et al.||1999| IDiemand fc Kuhlen||2008l |Springel et al.||2008l 



Stadel et al. 2009[ ), while observations indicate that only 



relatively few satellite galaxies exist around the Milky Way 
and M31 ( |Moore et al.|[l999} [Klypin et al.|[l999l ). Secondly, 
the highest resolution halo simulations show that the slope 
of the inner density profile decreases linearly at smaller radii 
(Navarro et al.| 1997| [Moore et al.|1999||Diemand et al.|2004 



Springel et al. 



2008 Stadel et al. 2009k , whereas the den- 



sity profiles inferred from galaxy rotation curves are signifi- 
cantly shallower ( Moore et al.|1999 l (and for recent studies 



see |Swaters et "a l 
|2008 | IGentile et 



2003 
al.||2009 



Salucci et aT^|2007l |de Blok et aT] 



and references there in). Thirdly, 
the observed number of dwarf galaxies in the voids appears 



to be far smaller than expected from CDM (Peebles 2001 



Tikhonov et al.||2009| [Peebles fc Nusser|[20"To[ ). Another ex 



ample is the excess in the prediction of dwarf galaxy con- 
centrations ( Lovell et al.|2011 l. Whilst, it has become clear 
that some of these discrepancies might be resolved through 
an improved understanding of galaxy formation, they have 
led some to consider changes to the ACDM paradigm. 

One possible solution might be warm dark matter 
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(WDM) (|Bond fc Szalay||1983| |Bardeen et al]|1986| |Bode 
|et al.|[2001| ). In this scenario, the dark particle is consid- 
ered to be lighter than its CDM counterpart, and so re- 
mains relativistic longer and also retains a thermal velocity. 
Since WDM particles are collisionless and decouple early, 
they may 'free-stream' or diffuse out of perturbations whose 
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size is smaller than the Jeans' lengtqjin the radiation dom- 
inated Universe (Kolb fc Turner|1 990[ ). This free-streaming 



of the WDM particles acts to damp structure formation on 
small scales. Two potential candidates are the sterile neu- 



trino ( Dodelson fc Widrow||l994 Shaposhnikov & Tkachev 
2006|), and the gravitino jEllis et al.|1984||Moroi et al.|1993 

tfanrooaW of al 1 0.0.71 In™-Kiimrvi7 o+ al lOflflSk KnfV. nf ™V,ioV 



Kawasaki et al.|[l997 Gorbunov et al.|2008 l, both of which 
require extensions of the standard model of particle physics. 

Recent observational constraints have suggested that 
sterile neutrinos can not be the dark matter: the Lyman 



alpha forest (Seljak et al. 2006 Boyarsky et al. 2009a) 



and QSO lensing (Miranda & Maccid 20071 bounds are 
m Vs > 8keV, whilst those from the X-ray background are 
m„. < 4keV ( |Boyarsky et al.||2008| p" 



However, a more re- 
cent assessment has suggested that a better motivated parti- 
cle physics model based on resonant production of the sterile 
neutrino, may evade these constraints: the Lyman alpha for- 
est bound is brought down to m Vs > 2keV and the X-ray 
background is pushed to m Vs < 50keV (for very low mixing 
angles) ( Boyarsky et al.|2009b l. It therefore seems that ad- 
ditional, independent methods for constraining the AWDM 
scenario would be valuable. 



In |Markovic et al.| ( |2010[ ) and |Smith fc Markovic| ( |201l[ ), 
it was proposed that the AWDM scenario could be tested 
through weak lensing by large-scale structure. The advan- 
tage of such a probe is that it is only sensitive to the total 
mass distribution projected along the line of sight. However, 
to obtain constraints on the WDM particle mass, an accu- 
rate model for the nonlinear matter clustering is required. 
In these papers, an approach based on the halo model was 
developed. Accurate predictions from this model require: de- 
tailed knowledge of the abundance of dark matter haloes, 
their spatial large-scale bias, and their density profiles. In 
these studies, it was assumed that the semi-analytic meth- 
ods, which were developed for CDM, would also apply to 
WDM. 

In this paper we perform a series of very high resolu- 
tion CDM and WDM A?-body simulations with the specific 
aim of exploring the halo model ingredients in the AWDM 
scenario. Over the past decade, there have been a limited 
number of numerical simulation studies of nonlinear struc- 



ture formation in the WDM model 


Colombi et al. 


1996 


Moore et al. 1999| |Coh'n et al.| 2000 


White & Croft 


2000 


Avila-Reese et al.|2001||Bode et al.|2001||Bullock et al. 
Zentner & Bullock|2003||Coh'n et al.|2008||Zavala et al. 


2002 
2009 


Maccid & Fontanot 


2010| |Lovell et al.|2011||Viel et al.;'2011 


Dunstan et al. 2011 


). In most of these previous studies, con- 



elusions have been drawn from object-by-object comparison 
of a relatively small number of haloes simulated in boxes of 
typical size L — 25 /i _1 Mpc. In this work we are more inter- 
ested in the overall impact that the WDM hypothesis has on 
the statistical properties of large-scale structures. We there- 
fore simulate boxes that are 10 times larger than have been 



1 Although originally defined in the context of gas dynamics, the 
Jeans length can be generalized to collisionless systems by replac- 
ing the sound speed with the velocity dispersion. The reason for 
this tight analogy lies in the linearized equation of perturbations, 
which has the same structure for gas and collisionless fluids (see 
|Peebles|1982| for more details). 

1 Lower bounds on the mass of a fully thermalized WDM particle 
can be obtained using Eq. ^ (see |Viel et al.|2005| . 



typically studied before, hence having roughly ~ 1000 times 
larger sampling volume. This means, that our conclusions 
will have greater statistical weight, than those from previous 
studies. Furthermore, our results should be less susceptible 
to finite volume effects, which can lead to underestimates of 
the nonlinear growth. 

The paper is structured as follows: In S|2]we provide a 
brief overview of the salient features of linear theory struc- 
ture formation in the WDM model and we review the halo 
model approach. In Sj3]we describe the A r -body simulations. 
In ij4] we explore the main ingredients of the halo model: 
the halo mass function, bias and density profiles. In i|5]we 
compare the halo model predictions for the matter power 
with our measurements from the simulations. Finally, in Sj6] 
we summarize our findings. 



2 THEORETICAL BACKGROUND 

In this section we summarize the linear theory for WDM 
and the nonlinear halo model in this framework. 



2.1 Linear theory evolution of WDM 

The physics of the free-streaming or diffusion of collision- 
less particles out o f dark matter pe r turba tions has been dis- 
cussed in detail by |Bond fc Szalay ( 1983[ P| An estimate for 
the free-streaming length can be obtained, by computing 
the comoving length scale that a particle may travel up un- 
til matter-radiation equality (£eq). At this point, the Jeans' 
length drops dramatically and perturbations may collapse 
under gravity. A simple formula for this is given by |Kolb"%t| 
Turnerl JT9901): 



v(t)dt 



cdt 



+ 



v(t)dt 



(1) 



where £nr is the epoch when the WDM particles become 
non-relativistic, which occurs when Twdm < towdmc 2 /3/cb, 
where Twdm and mwDM are the characteristic temperature 
and mass of the WDM particles. In the relativistic case, the 
mean peculiar velocity of the particle is simply v(t) ~ c. In 
the non-relativistic regime its momentum simply redshifts 
with the expansion: v oc a(t)~ . This leads to: 



1 , 1 1 * E Q 

1 H — log — — 

2 iNR 



(2) 



where th^nr) is the comoving size of the horizon at £nr- 
On inserting typical values for £nr we find the scaling: 



0.4 f m WDM \ 

V keV / 



-4/3 



liWDMB 

0.135 



[/i _1 Mpc](3) 



However, the real situation is more complex than this, since 
fluctuations inside the horizon grow logarithmically dur- 
ing radiation domination via the Meszaros effect and free- 
streaming does not switch off immediately after t^Q- To un- 
derstand the collisionless damping in more detail, one must 
numerically solve the coupled Einstein-Boltzmann system 
of equations for the various species of matter and radiation. 



3 For some recent theoretical treatments of WDM, also see |BoyJ, 
anovskyj ( [2010) 1 and |de Vega fc Sanchezj | |2010| |2uTT} 
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Several fitting formulae for the WDM density transfer func- 



tion have been proposed (Bardeen et al. 1986 Bode et al 



2001 1 and here we adopt the formula in Viel et al. (2005) 

1/2 



TwDivr(fc) = 



pCDM 

with (i — 1.12 as well as 
a = 0.049 



= [l + (afc) 2 "] 



-5/m 



(4) 



Mpc/h. (5) 



Note that in the above we are assuming that the WDM 



mwDM 


-l.n 


^WDM 


0.11 


■ h ■ 


L keV 




_ 0.25 _ 




.07. 



particle is fully thermalized. Following Viel et al. (20051, 
the masses of sterile neutrino WDM particles m Vs can be 
obtained from m-yvDM through the formula: 



m „ = 4.43keV 



/ m-WDM 

V IkeV 



4/3 



/ ffiWDM \ 

V 0.1225/ 



-1/3 



(6) 



The characteristic length-scale a is related to the free- 
streaming scale Af s , and we shall therefore make the def- 
inition that a = Af s ff is an effective free-streaming length 
scale. The length-scale Aff can be used to introduce the 
'free-streaming' mass scale: 



47T_ 



(7) 



where p is the background density of the universe. This mass 
scale is important as it defines the scale below which initial 
density perturbations are completely erased. 

We can define yet another length scale: the 'half-mode' 
length scale Ah m - This corresponds to the length scale at 
which the amplitude of the WDM transfer function is re- 
duced to 1/2. From Eq. Q we find: 

1/2 M 



A h 



2ttA 



eff fryn/S 



i 



13.93A?. 



(8) 



This length scale leads us to introduce another mass scale, 
the half-mode mass scale: 



Ah 



4tt_ / A 
-TP 



( Ahm \ ' 



2.7 x lO^Affs 



(9) 



This mass scale is where we expect the WDM to first affect 
the properties of dark matter haloes. In passing, this partly 



explains the claims made by Smith & Markovic ( 2011 1, that, 



for instance, the mass function of haloes would be signifi- 
cantly suppressed on mass scales M ~ 100Mf s . 

In Fig.[l]we show the relation between Mf B , Mh m and 
the mass of the WDM particle candidate for our adopted 
cosmological model. Three cases of relevance are apparent: 
M > Mhm, and haloes form hierarchically through accreting 
material; Mh m > M > Mf B and for these haloes the hierar- 
chy may fail with low mass haloes forming at the same time 
as higher mass haloes; finally Mf s > M no halo formation, 
unless through the fragmentation of larger structures. While 
the growth of overdensities is not affected above Mhm, it is 
suppressed between Mf s and Mhm, and should simply not 
take place below Mf s . 



2.2 Nonlinear evolution of WDM: the halo model 

Cosmological structure formation is a very complicated, 
highly nonlinear process that requires numerical simulation 
for a full understanding. However, the halo model approach 
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Figure 1. Free-streaming mass-scale (Mf s ) and half-mode mass 
scale (Mhm) as a function of the mass of the WDM particle 
(™wdm)- Haloes with masses M > Mh m) ma y form hierar- 
chically (upper right solid blue region). For haloes with masses 
Mh m > M > Mf s , hierarchical structure growth may fail (mid- 
dle green region). For haloes with masses M < Mf s , these may 
not form hierarchically since their initial peaks are completely 
erased (lower right empty region). However it is possible that 
such objects may emerge through fragmentation. The yellow dot- 
dash line denotes the current mwDM allowed by the Lyman al- 
pha forest (Boyarsky ct al. 2009a| (note that we have rescaled 
m Vs — > ?n wdm using Eq. (|6jl). 



gives a simplified analytical description of structure for- 
mation, which leads to surprisingly good results (Cooray 



& Shcth 2002, and references therein). Recently, the halo 
model has been adapted for the WDM cosmological model 



by Smith & Markovic (20111 and we now summarize their 
basic approach. 

The main idea of the halo model in WDM is to separate 
the density field into a halo component, adding up all bound 
structure, and a smooth component, standing for all matter, 
that has not collapsed due to free streaming. This is different 
to the standard approach of the CDM halo model, where all 
matter is supposed to be in bound structures. 

Thus the WDM density field has the form, 



p(x) = p a {x.) + ^2Miu(\x- 



,Mi 



(10) 



where p s is the smooth part of the density field and 
u(x, M) = ph(x\M)/M is the mass normalized density pro- 
file. The average densities of the smooth and the bound com- 
ponents are then given by 



(P) 



P = P s + Phi 



Ph = fP , 



(11) 



where / is the fraction of matter in bound objects. This can 
be calculated by integrating over the halo mass function 
weighted by halo mass: 



f== dlogMM-^— 



(12) 
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where dn — n(M)dM is the abundance of WDM haloes of 
mass M in the interval dM . The fraction / is equal to unity 
in a perfectly hierarchical universe and drops below unity 
as soon as the mass function is suppressed due to the free 
streaming. In a WDM universe the amount of suppression 
depends on the mass of the WDM particle. 

The power spectrum P(k) is defined by the relation 



<<5(k)J(k')> = (27r) 3 «S D (k + k')P(fc) 



(13) 



where So is the three dimensional Dirac delta function 
and 8(k) is the Fourier transform of the matter overden- 
sity <5(x) = (p(x) — p)/p- In terms of the different density 
components, we can write: 



5(k) = fh(k) + (1 - f)5 s (k), 



(14) 



where 5 X = (p x — p x )/p x with \ £ {h, s}. The power spec- 
trum of the halo model can now be determined by adding 
up the power spectra of the different density components as 
well as their cross terms, giving 

P(k) = (1 - f) 2 P BB (k) + 2(1 - /)/P sh (fc) + / 2 P hh (fc) . (15) 

The term Phh represents the power spectrum of matter 
trapped in haloes, the term P ss designates the power spec- 
trum of the smooth component and the term P a h denotes 
the cross-power spectrum between haloes and the smooth 
field. 

The term Phh can be separated into one- and two-halo 
terms, which describe the power coming from the same halo, 
and the one coming from distinct haloes, respectively. It can 
be expressed as: 



PhhW + fhhW 



Phh(fc) 

2 



dMi 



xP h c h (fe|Mi,M 2 



(16) 

(17) 
1 f°° 

Phh(fc) = ^2 dMn{M)M 2 u(k\M) , (18) 
Ph Jo 

where u(k\M) is the Fourier transform of the mass normal- 
ized density profile. In Eq. ( 17 1 we have introduced the power 
spectrum of halo centers P c {k\M\, M2), which in general 
is a complicated function of k and the halo masses Mi and 
M%. However, if we neglect halo exclusion and assume linear 
biasing with respect to the linear mass density, then we may 
write this as, 



P h c h (fc|M 1 ,M 2 ) ~ &i(Afi)&i(Af a )fli»(*) 



(19) 



In this case, the function is separable and this consider- 



ably simplifies the integrals in Eq. (17 1. This approximation 



breaks down on small, nonlinear scales, but on these scales, 
the two-halo term is sub-dominant. The error induced by 



this approximation ( 19 1 is most apparent at quasi-linear 
scales (k ~ [0.1, 1.0] /iMpc -1 ) and is < 30%. It is possible to 
lower this error to < 5% by using higher order perturbation 
theory techniques and by including halo exclusion (see for 
example |Smith et al.|[2011[ ). An easy but not fully consis- 
tent way of reducing the error down to < 10% is to do the 



following replacement in Equation (191: 

P lin (fc) -> P h aiofit(fc)W T H(fcP), R ~ 2/i~ 1 Mpc, 



(20) 
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Figure 2. Linear power spectra as a function of wavenumber 
in the CDM and WDM scenarios, at the initial redshift (z = 
49) of the simulations. Top panel: absolute dimensionless power: 
A 2 = k 3 P(k) /2-7T 2 . The lines denote the linear power spectrum 
where m^DM £ {00, 0.25, 0.5, 1.0} keV. Points denote the power 
spectra measured from the initial conditions of the N = 1024 3 
simulations. Bottom panel: ratio of the initial WDM and CDM 
power spectra. Points and lines unchanged. 



Phaioflt is the power spectrum calculated by the halofit 
code ( |Smith et aL|2003 l. 

The halo-smooth power spectrum is given by: 



Psh(fc) = 



Ph 



dMn(M)Mu(k\M)P£ s (k\M) 



(21) 



where P^ s (k\M) is the power spectrum of the halo centers 
with respect to the smooth mass field. On assuming that the 
smooth field and the halo density field are linearly biased 
with respect to the linear density field, we are lead to the 
relation: 



P s c h (fc|M) ~ b s b(M)P lin (k) 



(22) 



where b s is the linear bias of the smooth matter field defined 
in §4.3| Finally, the smooth field auto-power spectrum is 
given by 



P aa (k) = ^Plin(fc) 



(23) 



In order to reduce the error we can again replace the linear 
power spectrum in the Equations (|22[) and (|23[|, following 



the recipe of relation ( 20 1 



On combining these power spectra, weighted by the cor- 



rect functions of their mass fractions, a la Eq. ( 15 1, we find 



where Wth is the window function defined in §4.11 and 



the total halo model prediction for the nonlinear matter 
power spectrum in the WDM model. 
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Sim label 


m WDM [keV] 


Affc {h- L Mo} 




Llh^Mpc] 


A^part 


m p [h- 1 M @ } 
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256 3 


7.57 x 10 10 
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oo 








256 
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9.45 x 10 9 


10 


CDM-L 










1024 3 


1.18 x 10 9 


5 


WDM-1.25-S 










256 3 


7.57 x 10 iu 


20 


WDM-1.25-M 


1.25 


2.3 x 10 s 


6.3 x 10 9 
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512 3 


9.45 x 10 9 


10 


WDM-1.25-L 










1024 3 


1.18 x 10 9 


5 


WDM-l.O-S 










256 3 


7.57 x 10 10 


20 
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4.9 x 10 6 


1.3 x 10 10 
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9.45 x 10 9 
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1024 3 


1.18 x 10 9 
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WDM-0.75-S 
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7.57 x 10 1() 


20 


WDM-0.75-M 


0.75 


1.3 x 10 7 


3.4 x 10 10 


256 


512 3 


9.45 x 10 9 
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WDM-0.75-L 










1024 3 


1.18 x 10 9 


5 


WDM-0.5-S 










256 3 


7.57 x 10 10 


20 


WDM-0.5-M 


0.5 


4.9 x 10 7 


1.3 x 10 11 


256 


512 3 


9.45 x 10 9 


10 


WDM-0.5-L 










1024 3 


1.18 x 10 9 


5 


WDM-0.25-S 










256 3 


7.57 x 10 1() 


20 


WDM-0.25-M 


0.25 


5.0 x 10 8 


1.3 x 10 12 


256 


512 3 


9.45 x 10 9 


10 


WDM-0.25-L 










1024 3 


1.18 x 10 9 


5 



Table 1. WDM simulations. From left to right, columns are: simulation name (S=Small, M=Mcdium, L=Large); mass of WDM particle 
(ttiwdm); free-streaming mass-scale (Mf s ); half- mode mass-scale (Mh m ); simulation box-size (L); number of particles (iVpart); mass of 
simulation particles (m p ); comoving softening length (i so ft)- 



3 iV-BODY SIMULATIONS OF WDM 

In order to study nonlinear structure growth in the WDM 
model, we have generated a suite of iV-body simulations. 
These were executed on the zB0X3 supercomputer at the 
University of Zurich. Each simulation was performed us- 
ing PKDGRAV, a high order multipole tree-code with adaptive 
time stepping (Stadcl 20011. 

The cosmological parameters of the base ACDM model 
adopted, are consistent with the WMAP7 best-fit parame- 
ters ( jKomatsu et al.||2011[ ) and we take: the energy-density 
parameters in matter, vacuum energy and baryons to be 
O m = 0.2726, D.a = 0.7274, D. h = 0.046; the dimensionless 
Hubble parameter to be h = 0.704; the primordial power 
spectral index and present day normalization of fluctuations 
to be n s = 0.963, and a$ = 0.809. 

The CDM transfer function was generated using the 
code CAMB (Lewis et al.||2000l). The linear power spectrum 



for each WDM model was then obtained by multiplying the 
linear CDM power spectrum by T^ f - DM (k) from Eq. Ini- 
tial conditions for each WDM model were then generated 
at redshift z = 49 using the serial version of the publicly 
available 2LPT code ( Scoccimarro|1998 Crocce et al.|2006 |. 
In theory, we should also include a velocity dispersion due 
to the fact that the particles still retain a relic thermal ve- 
locity distribution. However, a quick calculation of the rms 
dispersion velocity, showed that these effects should be of 
marginal importance on scales > 50/i -1 kpc at the initial 
redshift, and of order > 1 /i -1 kpc at the present day for 
diwdm ^ 0.25 [keV]. We therefore assume that their inclu- 
sion will be a second order effect and so at this stage we 
neglect them. 

We generated initial conditions for a suite of simula- 
tions, one with a CDM particle and five with WDM particle 
masses mwDM G {0.25,0.5,0.75, 1.0, 1.25}keV. For all runs, 
we set the box length L = 256 h~ Mpc. This size is a com- 
promise between choosing a box small enough to accurately 
capture small-scale structure formation and large enough 
to confidently follow the linear evolution of the box-scale 



modes. This makes it possible for us to check agreement 
with the linear theory and to measure linear bias. 

Our simulations were also performed with three differ- 
ent mass resolutions: N = {256 3 , 512 3 , 1024 3 }. This enables 
us to differentiate between genuine structures and spurious 
structures, which can collapse out of the initial particle lat- 



tice (cf. |Wang fc White 2007 Polisensky fc Ricotti 20101 



Full details of the suite of simulations are summarized in 
Tabled] 

Dark matter haloes in the simulations were located us- 
ing the Friends-of- Friends algorithm ( Davis et al.|1985 I. We 
used a modified version of the skid code, with the linking 
length parameter set to the conventional value of b = 0.2. 

Fig. [2] compares the initial linear theory power spectra 
with the power spectra estimated from the initial conditions 
of the A^-body simulations, for the case N = 1024 3 . These 
results show that, at the initial time, the WDM linear theory 
distribution of fluctuations has been correctly seeded. It also 
shows a spike in the measured power spectrum at k — 8n 
which is a consequence of the initial particle distribution on 
a grid. 

Fig.[3]presents a pictorial view of the growth of structure 
in a selection of the simulations. The left column shows the 
density evolution in a slice through one of the CDM sim- 
ulations. The central and right panels show the same but 
for the case of WDM particles with biwdm = 1.0 keV and 
mwDM = 0.25 keV. From top to bottom the panels show 
results for z = 4.4, 1.1 and 0. 



4 HALO MODEL INGREDIENTS IN THE 
WDM SCENARIO 

In this section we detail the halo model ingredients and show 
how they change in the presence of our benchmark set of 
WDM particle masses. 
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Figure 3. Density maps from the N = 1024 simulations with a length of 50 h x Mpc and a depth of 2.5 h ^^Mpc. From top to bottom: 
z = 4.4, z = 1.1 and 2 = 0. From left to right: CDM, WDM with m p = l.OkeV and WDM with m p = 0.25keV. Whilst the WDM effects 
are barely discernible in the middle panels, they are very prominent in the right panels, where the voids are noticeably emptier than in 
CDM. 



4.1 Halo mass function 

In CDM the halo mass function can be explored through 



the excursion set formalism ( Press fc Schechter||1974 Bond 
et al.|1991| >: 



dn 1 p . .a! logo -2 

dlogM = ~2 M^dlogM 



(24) 



In the ellipsoidal collapse model of |Sheth fc Tormen| ( | 1999] ) , 
f(y) is given by 



fly) = M 



[l + (qu) 



-Pi g-9"/ 2 



s 2 c (t) 

a 2 {MY 



(25) 



with parameters: p = 0.3, q = 0.707 and with normalization 
parameter A — 0.3222. The linear theory collapse threshold 
is given by, 8 c (z) = 1.686/D(z), where D(z) is the linear 
theory growth function. The variance on mass scale M is, 



d j k 

(2^p 



Pun(k)W% H (kR) 



(26) 
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Figure 4. Measured mass function of the WDM simulations with m p = 0.25 keV and three different resolutions. The measurements lie 
below the Sheth-Tormen prediction, a well known result that is discussed in section]?] The upturn of the mass function due to artificial 
haloes is visible in the simulations of high and medium resolution. 



where Wth(v) = 3 [siny — y cos y] /y 3 and where the mass 
scale and radius of the filter function are related through 
the relation: M = 4nR 3 p/3. 

The main idea in the excursion set approach is that 
there is a monotonic mapping between the linear and non- 
linear density perturbations, averaged over a randomly se- 
lected patch of points in the space. Further, the mapping 
can be calculated using the spherical or ellipsoidal collapse 
approaches. The density perturbation in the patch will col- 
lapse to form a virialized object when the linearly extrapo- 
lated density in the patch reaches a certain collapse thresh- 
old. Despite the fact that this approach does not trace the 
full complexity of nonlinear structure formation, the actual 
predictions are in close agreement with measurements from 
simulations. That is, at least for a CDM cosmology with well 
calibrated values for the ellipsoidal parameters p and q and 
a given halo finding algorithm. One important assumption, 
which is implicit within this framework, is that structure 
formation must proceed hierarchically. 

In the WDM scenario, things are more complicated, 
since structure formation may not always proceed hierar- 
chically. As described in §j2j we can identify three regimes of 
interest: for M > Afhm, the variance of WDM fluctuations 
becomes virtually indistinguishable from that for CDM, and 
the excursion set approach should be valid; for M < Mf B 
all primordial overdensities are erased through diffusion of 
particles during the epoch of radiation domination and we 
expect that no hierarchical halo formation will take place 
on these mass scales. In between, where A/hm > M > Mf a , 
the WDM overdensity field is suppressed, but there is still 
some power left that may enable hierarchical collapse to take 



place. It is not clear a priori, how the mass function behaves 
on these scales and whether the extended Press-Schechter 
approach remains valid. We now investigate this using our 
simulations. 

In Fig.[4]we show the z = mass function of dark matter 
haloes for the case of rawDM = 0.25keV. The figure demon- 
strates the behaviour of the mass function as the simulation 
resolution is increased from N = 256 3 , to 512 3 , to 1024 3 
particles, denoted by the triangles, squares and circles, re- 
spectively. We can now see the effect of artificial clumping 
(cf. discussion in fj3J, which is manifest as the upturn of 
the curves at the low mass end of the mass function. One 
common approach to dealing with this artificial clumping is 
to assume that the simulations can be trusted down to the 
mass-scale just above the up-turn. We also find, in agree- 
ment with Wang & White (20071, that this mass-scale in- 
creases as N 1/J , i.e. the inter-particle spacing. In order to 
decrease the resolved mass by a factor of two, the particle 
resolution has to go up by a factor of eight. This is one of the 
main reasons why simulating WDM models is significantly 
more challenging than simulating CDM models. 

Fig. [4] also shows the prediction of the halo mass func- 
tion for CDM and for this WDM model, from the ST mass 
function. The figure clearly shows that the suppression of 
the ST model is not sufficiently strong to describe the data. 
In addition to this the ST mass function is diverging towards 
small masses, while we expect a realistic mass function to 
drop to zero at latest below the free streaming scale. 

Fig. [5] compares the measurements of the WDM mass 
functions from a selection of our highest resolution simula- 
tions with the CDM case. We note that, whilst for the case 
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Figure 5. Comparison between the Sheth-Tormen mass functions (black solid line for CDM, colored dashed lines for WDM) and 
the measurements from the simulations (black circles for CDM, colored squares, triangles and crosses for WDM). The grey solid lines 
correspond to the mass function fit of Eq. (271. 



of CDM the ST model is in very good agreement with the 
data, the WDM data all lie below the Sheth-Tormen predic- 
tion. That is, at least in the mass range above the artificial 
upturn of the mass function. 

Currently, there is no theoretical model that can ex- 
plain the discrepancy between the CDM and WDM mea- 
surements. We shall leave this as an issue for future study. 
However, it is possible to develop a fitting function that can 
describe the simulation results to high accuracy. As first 
noted in Smith & Markovic (2011 1, if one rescales the mass 
variable by Mf s , or equivalently by Mh m (cf. [J2J, then the 
mass functions for a wide variety of different values of mwDM 
all appear to fall upon the same locu^] 

In Fig.[6] we show that this scaling also works surpris- 
ingly well for the mass function measured from the simula- 
tions. We therefore look to fitting the rescaled mass func- 
tions. After trying various forms, we found that the function 



"WDM(Jtf) 

»wdm(^) 



= (1 + M hm /M)" 



(27) 



which has only one free parameter a = 0.6 was able to fit 
all of our data with an rms error well below five percent. 



4 We find that the locus of theory curves is much tighter than 
was first noted in |Smith fc Markovicj (2011| . This owes to the 
fact that, they adopted the free-streaming scale of |Bardeen et ah] 
(1986|;|Zentner fc Bull ock (2003), but used the transfer function 
of |Viel et al.| (2005| ) to generate the actual linear theory power 
spectra. This slight mismatch led to a slight off-set, which as Fig. [6] 
shows, is removed when consistent definitions for M[ s and 
are adopted. 
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Figure 6. Ratio between the WDM and CDM mass functions, 
as a function of halo mass, scaled in units of the half-mode mass- 
scale Mjjm. The measurements from the N = 1024 3 suite of WDM 
simulations are denoted by the point symbols. The dot-dashed 
lines denote the predictions from the |Sheth fc Tormen] (1999[ ) 
CDM mass function applied to WDM. The solid lines show the 
results from the fitting formula of Eq. (27|l. 
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Figure 7. Top panel: Comparison of halo bias in the CDM 
and mwDM = 0.25keV WDM model, as a function of wave num- 
ber. The open and solid points denote the results for CDM and 
WDM, respectively. Circles, stars, triangles and squares denote 
results for haloes with masses in the range: log 10 (M/[h _1 AfQ]) £ 
{[12.0,12.5], [12.5, 13.0], [13.0, 13.5] [> 13.5]} Lower panel: Ratio of 
the bias in the WDM model with that for the CDM model. 
For k > 0.1/iMpc - 1 we see a relative excess signal in the bias 
of haloes with M > 1O 12 /i _1 A/ in the WDM models. For 
k < 0.1/iMpc the trends are unclear owing to sample vari- 
ance. 



Note that in the above, n^ DM is the Sheth-Tormen model 
evaluated for the WDM model in question. The resulting 
mass functions are plotted as the grey solid lines in Fig. [5] 
A slightly worse fit may be obtained by using the function, 



4dm(M) 
"-cdm(M) 



(1 + M hm /M)" 



(28) 



with ft = 1.16 and this has the advantage that, one only 
needs to evaluate the CDM mass function and rescale the 
masses. We note that whilst this paper was being prepared, 



a similar study was presented by Dunstan et al. (2011 1, who 
showed that n^, DM provided a good description of their data 
but with the slightly higher value /3 = 1.2. 

Finally, we examined the evolution of the WDM mass- 



functions up to z — 1 and found that Eq. ( 27 1 also provides a 



good description of this data. The simplicity and generality 
of the fitting function Eq. \27\ is surprising and we think 
that it will be a useful empirical formula. 



4.2 Halo bias 

We are also interested in understanding how the density 
fields of dark matter haloes and matter are related in the 
WDM framework. This relation is usually termed bias, and 
if we assume that bias is local, deterministic, and linear, 
then we may write: 



5h(x|M) = 6(M)<5 m (x) , 



(29) 



where b(M) is the linear bias coefficient, which depends only 
on the mass of the halo. On using the excursion set formal- 
ism and the peak background split argument, one may ob- 
tain a prediction for b(M) (|Cole fc Kaiser|1989l|Mo fc White 
|1996||Sheth & Tormen||1999|): 



1 + 



qv 



2p 



S c (z) 8 c (z) [1 + (qv)P] 



(30) 



where the parameters p and q are as in Eq. (251. As was 



shown in Smith & Markovic (2011 1, if we apply this formula 



to the case of WDM, then we would expect to see that for 
M > Mh m the bias function is identical to that obtained for 
CDM. However, for M < Mh m we expect to find that the 
halo bias is increased relative to the CDM case. This occurs 
due to the fact that v tends towards a constant value for 
M < Mh m and so b ST becomes constant as well. We again 
use the simulations to investigate these predictions. 

In order to estimate the halo bias, we first sliced the 
halo distribution into a set of equal number density mass 
bins. Then, for each mass bin, we estimate the halo and 
matter auto-power spectra P hh (k\M) and P mm (fc), respec- 
tively. Our estimator for the bias at each fc-mode and in 
mass bin M a , can be written: 



bi(ki,M a ) 



P^(kj\M a ) - 

JD mm 



- l/n h (M a 



(31) 



where rih(M a ) is the number density of haloes for the mass 
bin a. 

Fig. [7] compares the scale-dependence of the halo bias, 
for several mass bins, and as a function of the wavemode, for 
the case of CDM (open points) and for the biwdm = 0.25 
keV WDM model (solid points). Note that here we actually 
present fe hm (fc) = p hm /p mm , where P hm is denned by the 
relation (27r) 3 5 D (k + k')P hm = (<5 h (k)<5 m (k')). In examining 
the ratio &wdm/&cdm, we see that there is increased bias in 
the WDM case. 

We then combine the estimates from each Fourier scale 
using a standard inverse variance weighted estimator (see 
e.g. Smith et al. 20071. Also, since, in this case, we are 



mainly interested in determining the effective linear bias, 
we only include modes with k < 0.1/iMpc -1 (cf. Fig. [7]). 
Fig.[8] presents the linear bias measurements together with 
the predictions from b ST (M) for a selection of the simulated 
WDM models. The four panels show the cases: CDM, top 
left; mwDM = l.OkeV, top right; jtiwdm = 0.5keV, bottom 
left; mwDM = 0.25keV, bottom right. 

Considering the high mass haloes, we find that the bias 
estimates for CDM and the WDM models appear to be in 
reasonable agreement with one another. At lower masses, 
however, there is a prominent increase in the bias for the 
WDM models with the lightest particle masses. We have 
found that, rather than a genuine effect due to WDM ini- 
tial conditions, this boost appears to be a manifestation of 
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Figure 8. Left panel: Linear halo bias for CDM (top left) and WDM (top right: 1.0 keV, bottom left: 0.5 keV and bottom right: 0.25 
keV). The filled and empty dots are measurements from the simulations with N = 1024 3 and TV = 512 3 , respectively. Error bars are 
calculated with an inverse variance wieghted estimator (see |Smith et al.|2007| . The solid lines correspond to the Sheth-Tormen model 
prediction of Eq. |30j. The linear halo bias of CDM is shown as a black dotted line for comparison. Right panel: Ratios between the 
WDM and the CDM linear halo bias for the N = 1024 3 runs. The error bars have been omitted for clarity. 



the artificial halo clumping discussed in ^ This becomes 
obvious by looking at the bottom-right plot in the panel 
on the left of Fig. [8] where the upturn of the high resolu- 
tion simulation (solid blue squares) is shifted with respect 
to the upturn in the low resolution simulation (empty blue 
squares). The mass scales of the upturn in the halo bias co- 
incides with the upturn in the mass function (see Fig.|4|. 
Importantly, this means that the halo-halo power spectrum 
is strongly contaminated by the spurious haloes, even on 
scales that are considered to be linear. To some extent this 
is not so surprising, given that below a certain mass scale 
we are dominated by spurious clumps. 

In general, when considering masses above the mini- 
mum mass-scale that we trust, rather than an excess bias 
with respect to 6 ST , we see that the estimates appear to lie 
slightly below the theoretical prediction at low masses. This 
has been observed in the CDM framework by |Tinker et al.| 
( |2010[ ) , and it seems to be the case for both our CDM and 
WDM simulations. However, for the case of mwDM = 0.25, 
we do note that, just above the non-physical upturn, there 
is a sign that bias in the WDM simulations is larger than in 
the CDM case. This trend is in qualitative agreement with 
the Sheth-Tormen prediction for WDM. However, the effect 
is small and of the order of the error bars and one would 
need both larger volume and higher resolution simulations 
to robustly confirm this. 



4.3 Anti-bias of the smooth component 

We also require the density field of the smooth matter. As 
for the case of the halo bias, if we assume that this is a 
linear, deterministic function of the matter density, then we 
may write the simple expression: 



S s (x) = b s 5(x) 



(32) 



where b B is the smooth bias parameter. As shown in |Smith] 
|fc Mar kovic ( 201l]) , this can be calculated using a mass con- 
servation argument, and one finds: 



/ 



- / dMMn(M)6i(M) 



< 1 



(33) 



Unlike the halo bias, which is mass dependent, the linear bias 
of the smooth component stays constant over all scales. In 
consequence, the smooth component of the power spectrum 
is directly proportional to the linear matter power spectrum. 
We shall leave it for future study to establish the veracity 
of this expression. 

4.4 Density profiles 

Over the years, extensive numerical work has shown that, 
for the case of the CDM model, the density profiles of dark 
matter haloes are reasonably well characterized by the NFW 



profile (Navarro et al. 1997 Moore et al. 
|et al.|2004||Springel et a,l.|2008[|Stadel et al 



1999 



Diemand 



2009). This has 



V 



(34) 



the universal form: 

P y(i + y) 2 

where the two parameters 5 B and r s represent a characteristic 
overdensity and scale radius. The mass of each halo can be 
determined by simply summing up the number of particles 
in a given object and multiplying by the particle mass. We 
can connect this to the virial radius through the relation: 

4-7T q 

M vir = — pAvirrvir, (35) 

where r v i r and A v i r are the virial radius and overdensity, 
respectively. The value of A v i r is typically chosen to de- 
note the overdensity for virialization, and here we adopt the 
value A vir = 200 (e.g. see |Sheth fc Tormen|1999[ ). However, 
the halo mass M v i r can also be obtained by integrating the 
density profile up to r v j r , which gives 



Mvir = <±-Kp5 s r a [log(l + C v ir) - Cvir/(1 + Cvir 



(36) 



where we have introduced the concentration parameter, de- 
fined as Cvir = r v i r /r B . On equating Eqs |35f and ( [36] ) we 
find that 
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Figure 9. Measurement of the halo profiles for CDM (black) and WDM (blue: 0.25 keV, green: 0.5 keV, red: 1.0 keV) for different halo 
masses. The profiles of each mass bin are coming from a randomly chosen halo, which is identified in the CDM and all WDM simulations. 



5 B = 



J A vir /3 



log(l + Cvir 



r/(l + C vir )] 



(37) 



This means that the original parameters {<5 s ,r s } of the 
NFW profile can be replaced by {M v i r ,c v ir}. Thus, given a 
simulated halo of mass M v i r , the model has one free param- 
eter, the concentration parameter c(M) 

Fig. [9] shows the density profiles of several randomly 
chosen haloes of different masses, for the case of CDM 
(black connected points). We have matched these objects 
with their counterpart haloes in the our standard set of 
WDM models and their profiles are also plotted (coloured 
connected points). While these profiles, on this logarithmic 
plot, all appear virtually indistinguishable for high masses, 
there does appear to be a net flattening off in the inner ra- 
dius for the galaxy mass haloes in the WDM model with 
mwDM = 0.25keV. One important point that can not be 
easily gleaned from the figure, is that there is an overall re- 
duction in the masses of all the smaller haloes. As we will 
see, this will have important consequences when we charac- 
terize the c(M) relation in the next section. 

Earlier work on this topic by |Moore et al.| ( |1999[ ) found 
that there was almost no perceptible difference between 
CDM haloes and haloes that formed from CDM initial con- 
ditions that had no small scale power below a certain scale. 



Subsequent work by Avila- Reese et al. ( 2001 1 and Colin 
et al. (20081, with a more careful treatment of the WDM 



transfer function, have shown more significant differences. 
However, in this case they were exploring models that were 
closer to HDM than WDM. We therefore conclude that our 
results are broadly consistent with all of these findings. 

One further point, is that for this small sample, we 
see no visible signs of the formation of a constant density 
core. This is in agreement with work Villacscusa-N avarro| 
& Dalai (20111. Adding thermal velocities into the simula- 



tions could in principle lead to the formation of a constant 
density core through the |Tremaine fc Gurrn] ( |1979[ ) limit on 
the fine grained phase space. However, the thermal veloci- 



ties of WDM cool down with the expansion of space and are 
already very small during the epoch of structure formation. 
Thus, if cores are induced, we expect that they will lie below 



the resolution limit of our simulations (see Kuzio de Naray 
et al.|2010| for more discussion of this) 



In summary, NFW profiles remain a valid approxima- 
tion for density profiles in our WDM simulations, given our 
spatial resolution and choice of mwDM- 

4.5 Concentration-mass relation 

As shown in the previous section the NFW model can be 
fully characterized by specifying the concentration mass re- 
lation. We now explore this for the case of WDM. 

In the CDM model, c v i r has been shown to be a mono- 



tonically decreasing function of M v i r (|Navarro et al. 


1997 


Bullock et al.|2001||Maccio et al.|2007[|2008||Neto et al. 


2007 


Prada et al.|2011 


1. One explanation for this, is that owing to 



the fact that haloes of different mass form at different times, 
they are therefore exposed to different background densities 
at collapse and this influences the final core overdensity. A 
denser background during collapse leads to generally higher 
concentrations. These ideas were encapsulated into a simple 
model for halo concentration by |Bullock et al.] ( |2001[ ): 

c vir = K(z c + l)/(z + 1) , (38) 

where z c is the redshift of collapse. This can be obtained by 
solving the relation a(M*,z) = 1.686, where M» = FM v i T , 
is defined to be a constant fraction of the virial mass. The 
two constants K and F must be calibrated using numerical 
simulations, and for our CDM simulations we found that 
K — 3.4 and F — 0.001 provided a good fit to the data. 
However, we note that the above arguments are only quali- 
tatively correct, since, as first pointed out by Bullock et al. 
(20011, there is a large scatter between c v i r and M v i r - This 



can, in part, be traced to the varying accretion histories and 
large-scale environments of different haloes of the same final 
mass. 
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Figure 10. Concentration to mass relation for CDM (top left) and for WDM with m = 1 keV (top right), m = 0.5 keV (bottom left) 
and m = 0.25 keV (bottom right). The colored symbols denote the median concentrations, while the dashed lines correspond to the 
Bullock model with F = 0.001 and K = 3.4. For comparison the Bullock model for CDM has been added to the WDM plots in form of 
a black dotted line. The gray lines are the 1-a contours. 



Turning to the WDM case, if we directly apply the Bul- 
lock model, but using the WDM linear power spectrum, then 
we find a suppression and a flattening of halo concentra- 
tions for masses M < Mhm- Similar to the mass function, 
this arises due to the fact that a(M) saturates to a constant 
value for masses approaching Mf s . We have tested these pre- 
dictions, by estimating the concentration parameters for all 
relaxed haloes in our CDM and WDM simulations that con- 
tain more than N — 500 particles (for full details of the 
method that we employ see |Macci6 et al.|2007| [2008] ). 

Fig. [To] shows the measured halo concentrations as a 
function of mass for a selection of the highest resolution 
CDM and WDM simulations. The gray solid lines corre- 
spond to the 1-a contours of the measurements, indicating a 
considerable spread in the concentration-mass relation. The 
large solid symbols denote the median, with the errors be- 
ing computed on the mean, i.e. we use a/y/Nl, where Ni 
are the number of haloes in the ith mass bin. The dashed 
lines denote the predictions from the Bullock model. For 
the CDM case it works reasonably well, especially with our 
modified parameters {K, F}. However, the model shows the 
wrong qualitative behavior for the WDM scenario: whilst 
the curve for the Bullock model always flattens out towards 
low masses, we see that for the cases of the lighter WDM 
particles, there is a turnover in the relation. This turnover 
in the c v i r -M relation at low masses is important, as it in- 
dicates the end of hierarchical collapse and the emergence 
of a period of top-down structure formation. As a test of 
these results we performed additional WDM runs with the 



cwdm(M) 
ccdm(M) 



= l+7i 



Mhn 

M 



(39) 



Gadget-2 gravity code (Springcl 2008) and we observe the 



same turnover in this independent set of simulations. 

In order to model the c v i r -M relation for WDM, we 
shall adapt the Bullock model. As in the case of the mass 
function, we do this by introducing a correction function 
described by the relation: 



where we have again rescaled the halo mass by M^ m . Least- 
squares optimization of the free-parameters gives: 71 = 15 
and 72 = 0.3. 

In Fig.JTTJ we compare the fitting function (gray solid 
lines) with the results from the simulations. The paramet- 
ric relation describes the cwum~M relation with a precision 
of better than 10% (the fit appears less satisfying for the 
case mwDM = 0.5keV, but only for the lower mass bins). 
Interestingly, the value 71 ~ 10, informs us that the c(M) 
relation is sensitive to the presence of WDM for mass scales 
one order of magnitude larger than for the mass function. 
As we will see in the next section, this will be important for 
modeling the nonlinear power spectrum on small scales. 



5 NONLINEAR POWER SPECTRUM 

5.1 Comparison with existing models 

In Fig. [12] we show the nonlinear matter power spectra esti- 
mated from our highest resolution CDM and WDM simula- 
tions. One can see that for k ^ 1 ft Mpc , there appears to 
be no obvious difference between the CDM and WDM mod- 
els under consideration. This is in stark contrast with the 
initial linear theory power spectra (cf. Fig.[2j), which show 
considerable damping for the same scales. Clearly nonlin- 
ear evolution has regenerated a high-fc tail to the power 
spectrum (cf. White fc Croft|20 00). At higher wavenumbers 
k > 1 h Mpc - , the situation is more interesting, and we see 
that the measured WDM power spectra are suppressed with 
respect to the CDM spectrum. The bottom panel of Fig. |13| 
quantifies this suppression in greater detail. Here we see that 
at k ~ 10 h Mpc -1 there is a 20% suppression in power for 
the case of mwDM = 0.25keV and this drops to ~ 2% for 
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Figure 11. Same as Fig. |10| but measurements are superimposed 
on one another and without error bars. The additional gray lines 
illustrate the fitting function from Eq. (|39| . 
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Figure 12. Nonlinear power spectra from the simulations (dots) 
and from the original halo model (lines), developed by |Smith fc| 
|Markovic| ( |2011} . Black corresponds to CDM and color to WDM 
(red: 1 keV, green: 0.5 keV, blue: 0.25 keV). The vertical gray 
dots indicate half the Nyquist frequency. 



the case mwDM = IkeV. The small difference between CDM 
and WDM at large scales (k < 1) is coming from a shift in 
the amplitude of the linear power spectrum, fixed with the 
same erg. 

We now explore whether the halo model approach, de- 
scribed in §2.2| can accurately reproduce our results for the 
WDM power spectra. In the original WDM halo model cal- 



culation of Smith & Markovic (20111, all of the model in- 



gredients (mass function, density profiles and halo bias re- 
lation), were obtained by assuming that the CDM relations 
also applied to the WDM case, provided one computes them 
using the appropriate linear power spectrum. The results of 
this approach are presented in Figs[l2]and[T3]as the coloured 
line styles. 

In Fi g.|13| we see that the halo model of | Smith fe| 
Markovic (2011 1 under-predicts the WDM power spectra by 
roughly ~ 10%. This is reasonably good, considering the as- 
sumptions that went into the model. This discrepancy was 
also noted in the study of Viel et al. (20111. In the bot- 
tom panel of Fig.[l3] we have also compared our nonlinear 
power spectr a with the predict ions from the fitting formula 
d201lb. For scales k < 10/iMpc" 1 



Viel et al. 



presented in 

we find that this fitting function provides an excellent de- 
scription of our data. However, for k > 10/iMpc -1 we find 
discrepancies, especially for the case mwDM = 0.25keV. 
Whether this is a genuine failing of the fitting formula is 
not clear, since this scale coincides with ~ kny/2, where 
&N y = nNgrid/L is the Nyquist frequency and we have used 
N gM = 2048. 

In summary, we find that the original halo model overes- 
timates the suppression of power due to WDM. This is not 
too surprising, since we have seen in the previous section 
that the original approximations for the halo mass function 
and concentrations turn out to be insufficient descriptions 
of the simulation data. 
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Figure 13. Top panel: Ratio of the simulated matter power spec- 
tra with respect to the halo model predictions as a function of 
wavenumber. Different coloured symbols denote the CDM and a 
selection of the WDM models. Bottom panel: Ratio of the WDM 
and CDM power spectra as a function of wavenumber. Points de- 
note the results from the ratios of simulation data; lines denote 
the halo model results. The gray solid lines correspond the fitting 
function from |Viel et al.| |2011| ). 



5.2 Towards an improved WDM halo model 

We now explore whether the halo model predictions can be 
improved by employing our better fitting functions for n(M) 
and p(r\M). Before making a final prediction for P(k), we 
first examine how each modification affects the predictions 
individually. 

If we implement our correction for the WDM mass func- 
tion in the P(k) predictions, then, since the abundance of 
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Figure 14. Ratio between modified versions of the halo model 
and the original version from |Smith fc Markovic] l |2011| |. The black 
solid line corresponds to CDM and the colored lines to WDM (red 
dotted: 1 keV, green dashed-dotted: 0.5 keV and blue dashed; 0.25 
keV). Top panel: only modification of the mass function. Middle 
panel: only modification of the concentration-mass relation. Bot- 
tom panel: modification of mass function and concentration-mass 
relation. 



small haloes is additionally suppressed with respect to the 
predictions of n ST for WDM, we should expect that there 
is an even stronger suppression in P(k). This conjecture is 
confirmed in the top panel of Fig. |14| which presents the 
ratio between the halo model with our modified mass func- 
tion and the original one. We clearly see that the ratio al- 
ways remains below unity. Somewhat surprisingly, we also 
note that a ~ 50% change in the abundance of 1O 12 /i _1 M0 
haloes, leads to a relatively small change, < 10%, in the 
power spectrum at k < 10/iMpc -1 . 

Next, if we instead implement our improved c v i T (M) 
relation, then we find that this has a more significant impact 
on the spectra. The central panel of Fig.[l4] shows the ratio 
between the halo model with the modified concentrations 
and the original one. We find that the suppression of the 
halo concentrations leads toa~ 50% boost in the power for 
k ~ 40hMpc _1 . 

The lower panel of Fig.[l4]shows the combined behavior 
of both corrections. The ratio between the fully modified 
halo model and the original one remains larger than unity. 
Thus, combination of the modified n(M) and c(M), leads to 
halo model predictions that have relatively more small-scale 
power. 

Finally, in Fig. [15] we present the comparison between 
our improved halo model and the nonlinear power spec- 
tra from the simulations. The top panel presents the ra- 
tios between the simulation data and the halo model predic- 
tions. The bottom panel shows the ratios of the WDM and 
CDM results for both the simulations and our modified halo 
model. Considering k > 3ftMpc -1 , whilst our modified halo 
model still has some problems predicting the overall abso- 
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Figure 15. Nonlinear power spectra from the simulations (dots) 
and from the fully corrected halo model (lines), including the fits 
for the mass function and the concentrations. The labeling is the 
same than in Fig. |12| The error of the halo model compared to the 
simulations has dropped below 10 percent (top pannel), the error 
on the ratio between WDM and CDM has dropped well below 5 
percent (bottom pannel). 



lute value of P(k), the relative changes between the WDM 
models and CDM are almost exactly predicted, they being 
accurate to better than ~ 2% down to k ~ 10/iMpc -1 . For 
scales beyond k > 10/iMpc -1 we see that the halo model 
also matches the simulations very well. However, again we 
note that these scales are beyond &N y /2 and so one might 
worry about aliasing effects. For k < 3 the error is below 
about 5%, this scales are however suffering from the diffi- 
culties in calculating P^ h , descirbed in j ]2.2| 

In summary, we conclude that our modified halo model 
is able to reproduce nonlinear WDM power spectra with the 
same accuracy as can currently be achieved for CDM. 



6 CONCLUSION 

In this paper we have explored nonlinear structure for- 
mation in the WDM cosmological model, through a large 
suite of cosmological iV-body simulations and through the 
halo model. The study was done for a set of fully ther- 
malized WDM models with particle masses in the set 
m — {0.25,0.5,0.75, 1.0, 1.25}keV. These masses range from 
purely pedagogical models, towards more realistic scenarios 
for the dark particle. 

For the simulations we chose a box size 
L = 256 7i _1 Mpc, which was small enough to resolve 
both the small scales, where WDM effects play an im- 
portant role, and the large scales, which are required for 
correct linear evolution of the box-modes. All models were 
simulated with N = {256 3 , 512 3 , 1024 3 } particles. This was 
done in order to disentangle physical effects from numerical 
ones. 

In the original halo model calculation for WDM by 



Smith & Markovic (20111, it was shown that in order to 
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make robust predictions, one requires good understanding of 
dark matter halo profiles, the mass function and halo bias. 
In this work we performed a detailed study of all of these 
ingredients. Our findings can be summarized as follows: 

(i) Mass function: Below a certain scale, the WDM mass 
function is suppressed with respect to CDM. This suppres- 
sion is considerably stronger than that obtained by simply 
applying the Sheth-Tormen approach together with the lin- 
ear power spectrum of WDM. In agreement with |Smith fe| 
Markovic (20111, we found that the mass functions for the 



quantities. Furthermore, we may also use the model to study 
the clustering of galaxies (Zehavi et al. 20051. It is hoped 



different WDM models could be transformed into a single 
locus of points. This was achieved by taking the ratio of the 
WDM mass function with that for CDM, and then rescaling 
the masses by Mh m (or equivalently Mf s ). We used a fitting 
function similar to that proposed in |Dunstan et al.] ( |2011[ | to 
link the Sheth-Tormen mass function to the measured one. 
The fitting function, which has only one free parameter, was 
able to reproduce all of the data with an accuracy of a few 
percent. We also found a strong boost in the mass function 
at very small mass scales. We showed that this was consis- 
tent with artificial halo formation around the initial particle 
lattice (cf. |Wang fc White|2007| . 

(ii) Halo bias: We measured the linear halo bias, using 
the four largest modes in our simulations. For smaller mass 
haloes, we found a small enhancement of the bias in WDM 
simulations, which was qualitatively consistent with the pre- 
dictions of |Smith fc Markovic|pOTT| . However, owing to the 
simulation box being too small, we were unable to quantify 
this more robustly. At very small masses we found a promi- 
nent boost in the bias. We found that this was again a sign 
of artificial halo formation. 

(iii) Density profiles: In the CDM model, the density pro- 
files of dark matter haloes can be characterized by an NFW 
profile, with a monotonically decreasing concentration-mass 
relation. In the WDM scenario, we have shown that the 
NFW profile remains valid for the models and resolution 
limits of our simulations, and we saw no evidence for 
a central density core. A simple adaption of the CDM 
concentration-mass relation, would suggest a strong flatten- 
ing towards small masses. Whilst, we found such a flatten- 
ing, the measurements in fact revealed a turnover towards 
smaller masses. This somewhat surprising result may be in- 
terpreted as a sign of top-down structure formation. We 
modelled the mean relation by adapting a fitting formula 
similar to that for the mass function. Our fit to the c(M) 
data was good to an accuracy of ~ 10%. Interestingly, we 
found that the deviations from CDM in the WDM model, 
appear in the c(M) relation for halo masses one order of 
magnitude larger than for the mass function. 

After analyzing these ingredients in detail and develop- 
ing new fitting functions for them, we were able to improve 
the small-scale performance of the WDM halo model. We 
found that for k > 3/iMpc" 1 , we could predict the absolute 
amplitude of the power spectrum to better than ~ 10%. 
However, we were able to predict the ratio of the WDM to 
CDM spectra, at better than < 2%. This was competitive 



with the latest fitting formulae ( Viel et al. 112011 1 



One of the many advantages of the halo model based 
approach, is that we may more confidently extrapolate our 
power spectra predictions to smaller scales than can be done 
from a fitting formula, since the model is built on physical 



that this may lead to a method for constraining WDM mod- 
els from galaxy clustering studies. Lastly, one further issue 
for future study, is to establish a better theoretical under- 
standing of what shapes the mass function and halo con- 
centrations in WDM. In particular, in finding the turnover 
in the concentration mass relation, have we really seen the 
reversal of bottom-up structure formation. This promises to 
be an interesting future challenge. 
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